//
// tetraTopology - class base for 3d tetrahedral topology
// By Antonio Miranda, Fev 2013
//////////////////////////////////////////////////////////////////////////

#include <math.h>
#include <algorithm>

#include "triface.h"
#include "../rtree/src/amr3bind.h"
#include "tetraTopology.h"


#define SPH_FACTOR                      05.0
#define SIGN(dot)                       (((dot)>0.0)?(1):(-1))
#define ABOUT_ZERO(val,tol)             ((val<tol)&&(val>(-tol)))


#define _NODE_BLOCK_ 100000


double tol_inters = 0.0;


//
//////////////////////////////////////////////////////////////////////////
tetraTopology::tetraTopology(void)
{
  m_nodeTree = NULL;
  m_faceTree = NULL;
  box.clear();

  // create trees
  m_nodeTree = RtreeCreate ( );
  m_faceTree = RtreeCreate ( );

  // 
  m_nodes.resize(_NODE_BLOCK_);

  m_smallestFace = -1.0;
  m_maxNodeId = -1;
}

//
//////////////////////////////////////////////////////////////////////////
tetraTopology::~tetraTopology(void)
{
  clear();
}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::clear()
{
  m_faces.clear();
  m_elems.clear();
  m_nodes.clear();


  m_faceMem.clear();
  m_elemMem.clear();
  m_nodeMem.clear();

  if (m_nodeTree != NULL)
    RtreeDestroy (m_nodeTree);
  if (m_faceTree != NULL)
    RtreeDestroy (m_faceTree);

  m_maxNodeId = -1;
}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::insertNode( int id, double x, double y, double z )
{
  node *new_node = m_nodeMem.get();
  new_node->coord[0] = x;
  new_node->coord[1] = y;
  new_node->coord[2] = z;
  new_node->real_id  = id;

  // update box
  box.update(x, y, z);

  // check capacity of nodes
  if (id > (int) m_nodes.capacity())
    m_nodes.resize(m_nodes.capacity()+_NODE_BLOCK_);

  // insert node
  m_nodes[id]= new_node;
  RtreeInsert (m_nodeTree, (void *) new_node, x, x, y, y, z, z);

  // update max id
  if (id > m_maxNodeId)
    m_maxNodeId = id;
}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::insertElement( int id, int n, int *conn, int gauss, int mat )
{
  int i, j;
  tetraelem *new_elem;
  double xmin, ymin, zmin, xmax, ymax, zmax;

  if (!(n == 4 || n == 10))
    return;

  // get new element
  new_elem = m_elemMem.get();

  // set vertex and id
  if (n == 4)
  {
    new_elem->vert[0] = conn[0] ;
    new_elem->vert[1] = conn[1] ;
    new_elem->vert[2] = conn[2] ;
    new_elem->vert[3] = conn[3] ;
    new_elem->type    = 0;
  }
  else
  {
    new_elem->vert[0] = conn[0] ;
    new_elem->vert[1] = conn[2] ;
    new_elem->vert[2] = conn[4] ;
    new_elem->vert[3] = conn[9] ;

    new_elem->vert[4] = conn[1] ;
    new_elem->vert[5] = conn[3] ;
    new_elem->vert[6] = conn[5] ;
    new_elem->vert[7] = conn[6] ;
    new_elem->vert[8] = conn[7] ;
    new_elem->vert[9] = conn[8] ;
    new_elem->type    = 1;
  }


  // insert new element
  new_elem->id = (int) m_elems.size();
  new_elem->gauss = gauss;
  new_elem->mat   = mat;
  m_elems.push_back(new_elem);

  // normal point to out of element
  int id_faces[4][3];
  new_elem->getFaces(id_faces);
  triface *faces[4] = {NULL, NULL, NULL, NULL};

  // try to find base 
  faces[0] = (triface *) m_nodes[new_elem->vert[0]]->faces.find(id_faces[0][0], id_faces[0][1], id_faces[0][2]);
  // try to find lateral faces
  for (i = 1; i < 4; ++i)
    faces[i] = (triface *) m_nodes[new_elem->vert[3]]->faces.find(id_faces[i][0], id_faces[i][1], id_faces[i][2]);

  // create new faces if necessary
  for (i = 0; i < 4; ++i)
  {
    if (faces[i] == NULL)
    {
      faces[i] = m_faceMem.get();
      faces[i]->addVertex(id_faces[i][0], id_faces[i][1], id_faces[i][2], m_nodes);

      // update adjacent faces to nodes
      for (j = 0; j < 3; ++j)
        m_nodes[id_faces[i][j]]->faces.insert (id_faces[i][0], id_faces[i][1], id_faces[i][2], 
        (void *) faces[i]);

      // insert face
      m_faces.push_back(faces[i]);
      faces[i]->getBox(xmin, ymin, zmin, xmax, ymax, zmax);
      RtreeInsert (m_faceTree, (void *)faces[i], xmin, xmax, ymin, ymax, zmin, zmax);

      // smallest faces
      double size = ((xmax-xmin) + (ymax-ymin) + (zmax-zmin));
      if ((m_smallestFace == -1.0 || size < m_smallestFace) && size != 0.0)
        m_smallestFace = size;
    }
  }

  // adjacent element <-> faces
  for (i = 0; i < 4; ++i)
  {
    // adjacent faces to element
    new_elem->faces[i] = faces[i];

    // adjacent element to faces
    if (faces[i]->adjElem[0] == -1)
      faces[i]->adjElem[0] = new_elem->id;
    else
      faces[i]->adjElem[1] = new_elem->id;

    // remove face from list if there is two element adjacent
    if (faces[i]->adjElem[1] != -1)
      faces[i]->boundary = false;
  }


}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::process()
{
  // boundary faces
  m_boundFaces.clear();
  size_t n = m_faces.size();
  for (size_t i = 0; i < n; ++i)
  {
    if (m_faces[i]->boundary)
      m_boundFaces.push_back(m_faces[i]);
  }

}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::printBoundFaces( char *file )
{

  FILE *ptr = NULL;
  
  fopen_s (&ptr, file, "wt");


  size_t n = m_boundFaces.size();

  fprintf (ptr,"Num = %d\n", n);

  for (size_t i = 0; i < n; ++i)
  {
    fprintf (ptr," %5d - %5d - %5d \n", m_boundFaces[i]->vert[0], m_boundFaces[i]->vert[1], m_boundFaces[i]->vert[2]);
  }
  
  fclose (ptr);
}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::insertNodalInf( string key, nodeInfo *info)
{
  m_nodalInf.insert(pair<string,nodeInfo *>(key,info));
}

//
//////////////////////////////////////////////////////////////////////////
triface * tetraTopology::getFace( int i, int j, int k )
{
  triface *face = NULL;

  // get node i
  node *node_i = m_nodes[i];
  face = (triface *) node_i->faces.find(i, j, k);

  return face;
}

//
//////////////////////////////////////////////////////////////////////////
void tetraTopology::getNodalInfo( string key, vector <nodeInfo *> &info)
{
  multimap<string,nodeInfo *>::iterator it;

  it = m_nodalInf.find(key);

  for (it = m_nodalInf.find(key); it != m_nodalInf.end(); ++it)
    info.push_back((*it).second);

}

